In [44]:
%reset
%pylab inline
In [45]:
import scipy.integrate as integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import math
from pycse import odelay
# Define stopping conditions, which can be used with odelay (from pycse)
def stop_yEquals0(state, t):
isterminal = True
direction = 0
value = state[1] # y = 0
return value, isterminal, direction
def stop_zEquals0(state, t):
isterminal = True
direction = 0
value = state[2] # z = 0
return value, isterminal, direction
# Not being used - don't have a decent test case with good values for the input state and BL1 (aka c2)
def linearDerivativesFunction(inputstate, timespan):
x, y, z, xdot, ydot, zdot = inputstate
#BL1 = 3.329168
#BL1 = 4.06107
BL1 = 0.012155092
derivs = [xdot,
ydot,
zdot,
2.0*ydot + (2.0*BL1 + 1.0)*x,
-2.0*xdot - (BL1 - 1)*y,
-BL1*z]
return derivs
# This full nonlinear version is what's being used
def nonlinearDerivativesFunction(inputstate, timespan):
x, y, z, xdot, ydot, zdot = inputstate
# distances
r1 = np.sqrt((mu+x)**2.0 + y**2.0 + z**2.0);
r2 = np.sqrt((1-mu-x)**2.0 + y**2.0 + z**2.0);
# masses
#m1 = 1 - mu;
#m2 = mu;
#G = 1;
derivs = [xdot,
ydot,
zdot,
x + 2.0*ydot + (1 - mu)*(-mu - x)/(r1**3.0) + mu*(1 - mu - x)/(r2**3.0),
y - 2.0*xdot - (1 - mu)*y/(r1**3.0) - mu*y/(r2**3.0),
-(1 - mu)*z/(r1**3.0) - mu*z/(r2**3.0)]
return derivs
In [46]:
#
#def generateClosedInitialConditions(inputstate, timespan, componentToVary, outputstate):
# propagate inputstate to plane crossing
#timespan, statesOverTime1, EventTime, EventState, EventIndex = odelay(nonlinearDerivativesFunction, inputstate, timespan, events=[stop_yEquals0])
#print 'initialstate1 = ', initialstate1
#print 't = ', t # timestamps corresponding to the output states over time
#print 'statesOverTime = ', statesOverTime
#print 'EventTime = ', EventTime
#print 'EventState = ', EventState
#print 'EventIndex = ', EventIndex
#print len(timespan)
#statesOverTime1 = integrate.odeint(nonlinearDerivativesFunction, initialstate1, timespan)
# if goal parameter is larger than tolerance, then
# adjust the componentToVary in the appropriate direction and repeat
In [47]:
# user inputs
ICset = 'Howell' # 'Sharp' 'Howell' 'Barbee'
ICtestcase = 0
numPoints = 2000
# create a dictionary for the initial conditions
ICs = dict()
# From Sharp, A Collection of Restricted Three-Body Test Problems
# For problems 1 to 15, mu = 0.012277471 and for problems 16 to 20, mu = 0.000953875
ICs['Sharp'] = {'mu': np.ones(20),
'X': np.zeros(20),
'Z': np.zeros(20),
'Ydot': np.zeros(20),
'T': np.zeros(20)}
Sharp_X_Z_Ydot_T = np.matrix([[0.994000E+00, 0.0, -0.21138987966945026683E+01, 0.54367954392601899690E+01],
[0.994000E+00, 0.0, -0.20317326295573368357E+01, 0.11124340337266085135E+02],
[0.994000E+00, 0.0, -0.20015851063790825224E+01, 0.17065216560157962559E+02],
[0.997000E+00, 0.0, -0.16251217072210773125E+01, 0.22929723423442969481E+02],
[0.879962E+00, 0.0, -0.66647197988564140807E+00, 0.63006757422352314657E+01],
[0.879962E+00, 0.0, -0.43965281709207999128E+00, 0.12729711861022426544E+02],
[0.879962E+00, 0.0, -0.38089067106386964470E+00, 0.19138746281183026809E+02],
[0.997000E+00, 0.0, -0.18445010489730401177E+01, 0.12353901248612092736E+02],
[0.100000E+01, 0.0, -0.16018768253456252603E+01, 0.12294387796695023304E+02],
[0.100300E+01, 0.0, -0.14465123738451062297E+01, 0.12267904265603897140E+02],
[0.120000E+01, 0.0, -0.71407169828407848921E+00, 0.18337451820715063383E+02],
[0.120000E+01, 0.0, -0.67985320356540547720E+00, 0.30753758552146029263E+02],
[0.120000E+01, 0.0, -0.67153130632829144331E+00, 0.43214375227857454128E+02],
[0.120000E+01, 0.0, -0.66998291305226832207E+00, 0.55672334134347612727E+02],
[0.120000E+01, 0.0, -0.66975741517271092087E+00, 0.68127906604713772763E+02],
[-0.102745E+01, 0.0, 0.40334488290490413053E-01, 0.18371316400018903965E+03],
[-0.976680E+00, 0.0, -0.61191623926410837000E-01, 0.17733241131524483004E+03],
[-0.766650E+00, 0.0, -0.51230158665978820282E+00, 0.17660722897242937108E+03],
[-0.109137E+01, 0.0, 0.14301959822238380020E+00, 0.82949461922342093092E+02],
[-0.110137E+01, 0.0, 0.15354250908611454510E+00, 0.60952121909407746612E+02]])
ICs['Sharp']['mu'][0:15] *= 0.012277471
ICs['Sharp']['mu'][15:20] *= 0.000953875
ICs['Sharp']['X'] = np.array(Sharp_X_Z_Ydot_T[:,0])
ICs['Sharp']['Z'] = np.array(Sharp_X_Z_Ydot_T[:,1])
ICs['Sharp']['Ydot'] = np.array(Sharp_X_Z_Ydot_T[:,2])
ICs['Sharp']['T'] = np.array(Sharp_X_Z_Ydot_T[:,3])
# From Howell, Three-Dimensional, Periodic, 'Halo' Orbits
ICs['Howell'] = {'mu': [0.04, 0.04],
'X': [0.723268, 0.723268],
'Z': [0.040000, -0.040000],
'Ydot': [0.198019, 0.198019],
'T': [1.300177*2.0, 1.300177*2.0]}
# From Barbee, Notional Mission 4 (Earth-Moon)
ICs['Barbee'] = {'mu': [0.012277471],
'X': [0.862307159058101],
'Z': [0.0],
'Ydot': [-0.187079489569182],
'T': [2.79101343456226]}
# assign simulation variables using specified elements from dictionary of IC's
mu = ICs[ICset]['mu'][ICtestcase]
timespan = np.linspace(0, ICs[ICset]['T'][ICtestcase], numPoints)
initialstate1 = [ICs[ICset]['X'][ICtestcase], 0, ICs[ICset]['Z'][ICtestcase],
0, ICs[ICset]['Ydot'][ICtestcase], 0]
In [48]:
# Compute the locations of the libration points
# ported this from Matlab code found at http://www.math.rutgers.edu/~jmireles/matLabPage.html - need to QA against a textbook
# In numpy.roots(p), p[0] corresponds to highest-order term in polynomial
l = 1.0 - mu
# L1
p_L1 = [1, 2.0*(mu-l), l**2.0 - 4.0*l*mu + mu**2.0, 2.0*mu*l*(l-mu) + mu-l, mu**2.0*l**2.0 + 2.0*(l**2.0 + mu*2.0), mu**3.0 - l**3.0]
L1 = np.real([i for i in np.roots(p_L1) if (i > -mu and i < l)])
L1 = np.append(L1, [0.0, 0.0])
# L2
p_L2 = [1, 2.0*(mu-l), l**2.0 - 4.0*l*mu + mu**2.0, 2*mu*l*(l-mu) - (mu+l), mu**2.0*l**2.0 + 2.0*(l**2.0 - mu**2.0), -(mu**3.0 + l**3.0)]
L2 = np.real([i for i in np.roots(p_L2) if (i > -mu and i > l)])
L2 = np.append(L2, [0.0, 0.0])
# L3
p_L3 = [1.0, 2.0*(mu-l), l**2.0 - 4.0*mu*l + mu**2.0, 2.0*mu*l*(l-mu) + (l+mu), mu**2.0*l**2.0 + 2*(mu**2 - l**2.0), l**3.0 + mu**3.0]
L3 = np.real([i for i in np.roots(p_L3) if i < -mu])
L3 = np.append(L3, [0.0, 0.0])
# L4
L4 = [-mu + 0.5, np.sqrt(3.0)/2.0, 0.0]
# L5
L5 = [-mu + 0.5, -np.sqrt(3.0)/2.0, 0.0]
print L1
print L2
print L3
print L4
print L5
In [49]:
# integrate first satellite
statesOverTime1 = integrate.odeint(nonlinearDerivativesFunction, initialstate1, timespan)
#timespan, statesOverTime1, EventTime, EventState, EventIndex = odelay(nonlinearDerivativesFunction, initialstate1, timespan, events=[stop_zEquals0])
#print 'initialstate1 = ', initialstate1
#print 't = ', t # timestamps corresponding to the output states over time
#print 'statesOverTime = ', statesOverTime
#print 'EventTime = ', EventTime
#print 'EventState = ', EventState
#print 'EventIndex = ', EventIndex
#print len(timespan)
x1, y1, z1, xdot1, ydot1, zdot1 = statesOverTime1.T
# identify x-coordinate corresponding to maximum y-amplitude,
# y-coordinate corresponding to maximum x-amplitude,
# and z-coordinate corresponding to maximum y-amplitude #y=0
center = [x1[y1.argmax()], y1[x1.argmax()], z1[y1.argmax()]] #np.abs(y1).argmin()]]
print 'center = ', center
In [50]:
# plot in RLP
figRLP, ((axXZ, axYZ), (axXY, ax3D)) = plt.subplots(2, 2)
figRLP.suptitle('Offset between Satellites 1 and 2 in RIC Frame')
figRLP.set_size_inches(10, 10)
figRLP.subplots_adjust(hspace=0.2)
figRLP.subplots_adjust(wspace=0.5)
# XZ Plane
axXZ.plot(x1, z1, 'r-')
axXZ.plot(L1[0], L1[2], 'g*', markersize=10)
axXZ.plot(center[0], center[2], 'w*', markersize=10)
axXZ.set_title('XZ Plane')
axXZ.xaxis.set_label_text('X axis')
axXZ.yaxis.set_label_text('Z axis')
axXZ.set_aspect('equal')
# YZ Plane
axYZ.plot(y1, z1, 'r-')
axYZ.plot(L1[1], L1[2], 'g*', markersize=10)
axYZ.plot(center[1], center[2], 'w*', markersize=10)
axYZ.set_title('YZ Plane')
axYZ.xaxis.set_label_text('Y axis')
axYZ.yaxis.set_label_text('Z axis')
axYZ.set_aspect('equal')
# XY Plane
axXY.plot(x1, y1, 'r-')
axXY.plot(L1[0], L1[1], 'g*', markersize=10)
axXY.plot(center[0], center[1], 'w*', markersize=10)
axXY.set_title('XY Plane')
axXY.xaxis.set_label_text('X axis')
axXY.yaxis.set_label_text('Y axis')
axXY.set_aspect('equal')
# 3D View
ax3D.axis('off')
ax3D = figRLP.add_subplot(224, projection='3d')
ax3D.plot(x1, y1, z1, 'ro-')
ax3D.plot([L1[0]], [L1[1]], [L1[2]], 'g*', markersize=10)
ax3D.plot([center[0]], [center[1]], [center[2]], 'w*', markersize=10)
ax3D.set_title('3D View in RLP Frame')
ax3D.xaxis.set_label_text('X axis')
ax3D.yaxis.set_label_text('Y axis')
ax3D.zaxis.set_label_text('Z axis')
#ax3D.set_aspect('equal')
#figXY, axXDotYDot = plt.subplots()
#axXDotYDot.plot(timespan, xdot)
#axXDotYDot.plot(timespan, ydot)
Out[50]:
In [51]:
# build RIC frame based on satellite 1
rVec = np.array([x1-center[0], y1-center[1], z1-center[2]]).T
vVec = np.array([xdot1, ydot1, zdot1]).T
cVec = np.cross(rVec, vVec)
iVec = np.cross(cVec, rVec)
# unitize RIC frame vectors
for ii in range(0, len(rVec)):
rVec[ii] /= np.linalg.norm(rVec[ii])
vVec[ii] /= np.linalg.norm(vVec[ii])
cVec[ii] /= np.linalg.norm(cVec[ii])
iVec[ii] /= np.linalg.norm(iVec[ii])
# build VNB frame based on satellite 1
#rVec = np.array([x1-center[0], y1-center[1], z1-center[2]]).T
#vVec = np.array([xdot1, ydot1, zdot1]).T
nVec = np.cross(rVec, vVec)
bVec = np.cross(vVec, nVec)
# unitize VNB frame vectors
for ii in range(0, len(rVec)):
rVec[ii] /= np.linalg.norm(rVec[ii])
vVec[ii] /= np.linalg.norm(vVec[ii])
nVec[ii] /= np.linalg.norm(nVec[ii])
bVec[ii] /= np.linalg.norm(bVec[ii])
In [52]:
## Theta Angle Approach
# compute angle theta (in XY plane for now) wrt center for each point along the orbit
theta = np.zeros(len(x1))
for ii in range(0, len(x1)):
theta[ii] = math.degrees( math.atan2( x1[ii] - center[0], y1[ii] - center[1] ))
# user input:
thetadirection = np.sign(theta[1] - theta[0])
desiredtheta = theta[50]# + thetadirection*5 # theta[2] # degrees
figTheta, axTheta = plt.subplots()
axTheta.plot(theta, 'o-')
figTheta.suptitle('Theta over time of Satellite 1')
# copy state from desired 'theta' value into second satellite initial state
desiredthetaindex = (np.abs(theta-desiredtheta)).argmin()
print 'theta[0] = ', theta[0]
print 'index = ', desiredthetaindex
print 'theta[index] = ', theta[desiredthetaindex]
initialstate2 = [x1[desiredthetaindex], y1[desiredthetaindex], z1[desiredthetaindex],
xdot1[desiredthetaindex], ydot1[desiredthetaindex], zdot1[desiredthetaindex]]
## VNB velocity offset approach
# compute a small offset to velocity along the velocity direction
# assign initial state to second satellite
#initialstate2 = [x1[0], y1[0], z1[0],
# xdot1[0]*1.001, ydot1[0]*1.001, zdot1[0]*1.001]
#generate closed halo...
## VNB position offset approach
# compute a small offset to position along the velocity direction
# assign initial state to second satellite
#initialstate2 = [x1[0] + xdot1[0]*0.001, y1[0] + ydot1[0]*0.001, z1[0] + zdot1[0]*0.001,
# xdot1[0], ydot1[0], zdot1[0]]
In [53]:
# integrate second satellite with same timespan as first satellite
statesOverTime2 = integrate.odeint(nonlinearDerivativesFunction, initialstate2, timespan)
x2, y2, z2, xdot2, ydot2, zdot2 = statesOverTime2.T
# compute trajectory offset in RLP frame
dx = x1 - x2
dy = y1 - y2
dz = z1 - z2
dxdot = xdot1 - xdot2
dydot = ydot1 - ydot2
dzdot = zdot1 - zdot2
# compute total distance offset
distance = np.zeros(len(dx))
for ii in range(0, len(dx)):
distance[ii] = np.linalg.norm([dx[ii], dy[ii], dz[ii]])
# plot total distance offset over time
figDeltaMag, axDeltaMag = plt.subplots()
axDeltaMag.plot(timespan, distance, 'o-')
figDeltaMag.suptitle('Total Distance Offset Over Time')
#figXY, axXDotYDot = plt.subplots()
#axXDotYDot.plot(timespan, xdot)
#axXDotYDot.plot(timespan, ydot)
Out[53]:
In [54]:
# plot both satellites in RLP frame
figRLP2, ((axXZ2, axYZ2), (axXY2, ax3D2)) = plt.subplots(2, 2)
figRLP2.suptitle('Satellites 1 and 2 in RLP Frame')
figRLP2.set_size_inches(10, 10)
figRLP2.subplots_adjust(hspace=0.2)
figRLP2.subplots_adjust(wspace=0.5)
# XZ Plane
axXZ2.plot(x1, z1, 'r-')
axXZ2.plot(x2, z2, 'b-')
axXZ2.plot(L1[0], L1[2], 'g*', markersize=10)
axXZ2.plot(center[0], center[2], 'w*', markersize=10)
axXZ2.set_title('XZ Plane')
axXZ2.xaxis.set_label_text('X axis')
axXZ2.yaxis.set_label_text('Z axis')
axXZ2.set_aspect('equal')
# YZ Plane
axYZ2.plot(y1, z1, 'r-')
axYZ2.plot(y2, z2, 'b-')
axYZ2.plot(L1[1], L1[2], 'g*', markersize=10)
axYZ2.plot(center[1], center[2], 'w*', markersize=10)
axYZ2.set_title('YZ Plane')
axYZ2.xaxis.set_label_text('Y axis')
axYZ2.yaxis.set_label_text('Z axis')
axYZ2.set_aspect('equal')
# XY Plane
axXY2.plot(x1, y1, 'r-')
axXY2.plot(x2, y2, 'b-')
axXY2.plot(L1[0], L1[1], 'g*', markersize=10)
axXY2.plot(center[0], center[1], 'w*', markersize=10)
axXY2.set_title('XY Plane')
axXY2.xaxis.set_label_text('X axis')
axXY2.yaxis.set_label_text('Y axis')
axXY2.set_aspect('equal')
# 3D view
ax3D2.axis('off')
ax3D2 = figRLP2.add_subplot(224, projection='3d')
ax3D2.plot(x1, y1, z1, 'ro-')
ax3D2.plot(x2, y2, z2, 'bo-')
ax3D2.plot([L1[0]], [L1[1]], [L1[2]], 'g*', markersize=10)
ax3D2.plot([center[0]], [center[1]], [center[2]], 'w*', markersize=10)
ax3D2.set_title('3D View in RLP Frame')
ax3D2.xaxis.set_label_text('X axis')
ax3D.yaxis.set_label_text('Y axis')
ax3D.zaxis.set_label_text('Z axis')
#ax3D2.set_aspect('equal')
Out[54]:
In [55]:
# plot offset in RLP frame
figDeltaRLP, ((axDeltaXZ, axDeltaYZ), (axDeltaXY, ax3DDelta)) = plt.subplots(2, 2)
figDeltaRLP.suptitle('Offset between Satellites 1 and 2 in RLP Frame')
figDeltaRLP.set_size_inches(10, 10)
figDeltaRLP.subplots_adjust(hspace=0.2)
figDeltaRLP.subplots_adjust(wspace=0.5)
# XZ Plane
axDeltaXZ.plot(dx, dz, 'b-')
axDeltaXZ.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaXZ.set_title('XZ Plane')
axDeltaXZ.xaxis.set_label_text('X axis')
axDeltaXZ.yaxis.set_label_text('Z axis')
#axDeltaXZ.set_aspect('equal')
# YZ Plane
axDeltaYZ.plot(dy, dz, 'b-')
axDeltaYZ.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaYZ.set_title('YZ Plane')
axDeltaYZ.xaxis.set_label_text('Y axis')
axDeltaYZ.yaxis.set_label_text('Z axis')
#axDeltaYZ.set_aspect('equal')
# XY Plane
axDeltaXY.plot(dx, dy, 'b-')
axDeltaXY.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaXY.set_title('XY Plane')
axDeltaXY.xaxis.set_label_text('X axis')
axDeltaXY.yaxis.set_label_text('Y axis')
#axDeltaXY.set_aspect('equal')
# plot offset in RLP frame
ax3DDelta.axis('off')
ax3DDelta = figDeltaRLP.add_subplot(224, projection='3d') # "224" means "2x2 grid, 4th subplot"
ax3DDelta.plot(dx, dy, dz, 'b-')
ax3DDelta.plot([0.0], [0.0], [0.0], 'ro', markersize=10)
ax3DDelta.set_title('3D View of Offset in RLP Frame')
ax3DDelta.xaxis.set_label_text('X axis')
ax3DDelta.yaxis.set_label_text('Y axis')
ax3DDelta.zaxis.set_label_text('Z axis')
#ax3DDelta.set_aspect('equal')
Out[55]:
In [56]:
# compute trajectory offset in RIC frame
dr = np.zeros(len(dx))
di = np.zeros(len(dx))
dc = np.zeros(len(dx))
for ii in range(0, len(rVec)):
dr[ii] = np.dot([dx[ii], dy[ii], dz[ii]], rVec[ii])
di[ii] = np.dot([dx[ii], dy[ii], dz[ii]], iVec[ii])
dc[ii] = np.dot([dx[ii], dy[ii], dz[ii]], cVec[ii])
print [dx[ii], dy[ii], dz[ii]]
print rVec[ii]
print np.dot([dx[ii], dy[ii], dz[ii]], rVec[ii])
print dr[ii], di[ii], dc[ii]
print len(dr), len(di)
# plot offset in RIC frame
figDeltaRIC, ((axDeltaRC, axDeltaIC), (axDeltaRI, ax3DRIC)) = plt.subplots(2, 2)
figDeltaRIC.suptitle('Offset between Satellites 1 and 2 in RIC Frame')
figDeltaRIC.set_size_inches(10, 10)
figDeltaRIC.subplots_adjust(hspace=0.2)
figDeltaRIC.subplots_adjust(wspace=0.5)
# RC Plane
axDeltaRC.plot(dr, dc, 'b-')
axDeltaRC.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaRC.set_title('RC Plane')
axDeltaRC.xaxis.set_label_text('R axis')
axDeltaRC.yaxis.set_label_text('C axis')
#axDeltaRC.set_aspect('equal')
# IC Plane
axDeltaIC.plot(di, dc, 'b-')
axDeltaIC.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaIC.set_title('IC Plane')
axDeltaIC.xaxis.set_label_text('I axis')
axDeltaIC.yaxis.set_label_text('C axis')
#axDeltaIC.set_aspect('equal')
# RI Plane
axDeltaRI.plot(dr, di, 'b-')
axDeltaRI.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaRI.set_title('RI Plane')
axDeltaRI.xaxis.set_label_text('R axis')
axDeltaRI.yaxis.set_label_text('I axis')
#axDeltaRI.set_aspect('equal')
# plot offset in RIC frame
ax3DRIC.axis('off')
ax3DRIC = figDeltaRIC.add_subplot(224, projection='3d') # "224" means "2x2 grid, 4th subplot".
ax3DRIC.plot(dr, di, dc, 'b-')
ax3DRIC.plot([0.0], [0.0], [0.0], 'ro', markersize=10)
ax3DRIC.set_title('3D View of Offset in RIC Frame')
ax3DRIC.xaxis.set_label_text('R axis')
ax3DRIC.yaxis.set_label_text('I axis')
ax3DRIC.zaxis.set_label_text('C axis')
#ax3DRIC.set_aspect('equal')
Out[56]:
In [57]:
# compute trajectory offset in VNB frame
dv = np.zeros(len(dx))
dn = np.zeros(len(dx))
db = np.zeros(len(dx))
for ii in range(0, len(rVec)):
dv[ii] = np.dot([dx[ii], dy[ii], dz[ii]], vVec[ii])
dn[ii] = np.dot([dx[ii], dy[ii], dz[ii]], nVec[ii])
db[ii] = np.dot([dx[ii], dy[ii], dz[ii]], bVec[ii])
#plt.xkcd(scale=1, length=100, randomness=2) # fun
#plt.rcdefaults()
# plot offset in VNB frame
figDeltaVNB, ((axDeltaVB, axDeltaNB), (axDeltaVN, ax3DVNB)) = plt.subplots(2, 2)
figDeltaVNB.suptitle('Offset between Satellites 1 and 2 in VNB Frame')
figDeltaVNB.set_size_inches(10, 10)
figDeltaVNB.subplots_adjust(hspace=0.2)
figDeltaVNB.subplots_adjust(wspace=0.5)
# VB Plane
axDeltaVB.plot(dv, db, 'b-')
axDeltaVB.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaVB.set_title('VB Plane')
axDeltaVB.xaxis.set_label_text('V axis')
axDeltaVB.yaxis.set_label_text('B axis')
#axDeltaVB.set_aspect('equal')
# NB Plane
axDeltaNB.plot(dn, db, 'b-')
axDeltaNB.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaNB.set_title('NB Plane')
axDeltaNB.xaxis.set_label_text('N axis')
axDeltaNB.yaxis.set_label_text('B axis')
#axDeltaNB.set_aspect('equal')
# VN Plane
axDeltaVN.plot(dv, dn, 'b-')
axDeltaVN.plot(0.0, 0.0, 'ro', markersize=10)
axDeltaVN.set_title('VN Plane')
axDeltaVN.xaxis.set_label_text('V axis')
axDeltaVN.yaxis.set_label_text('N axis')
#axDeltaVN.set_aspect('equal')
# plot offset in VNB frame
ax3DVNB.axis('off')
ax3DVNB = figDeltaVNB.add_subplot(224, projection='3d') # "224" means "2x2 grid, 4th subplot"
ax3DVNB.plot(dv, dn, db, 'b-')
ax3DVNB.plot([0.0], [0.0], [0.0], 'ro', markersize=10)
ax3DVNB.set_title('3D View of Offset in VNB Frame')
ax3DVNB.xaxis.set_label_text('V axis')
ax3DVNB.yaxis.set_label_text('N axis')
ax3DVNB.zaxis.set_label_text('B axis')
#ax3D.set_aspect('equal')
Out[57]:
In [57]: